Research Question: Do damselfish (Dascyllus flavicaudus) enhance coral skeletal growth rates in Pocillopora spp., and does this effect vary with wound size?
Experimental Design:
This document contains:
Key Methodology: We use allometric scaling to correct for initial coral size differences, allowing fair comparison of growth rates across all fragments.
# ==================== Core Packages ====================
library(tidyverse) # Data wrangling and visualization
library(here) # Relative file paths
library(janitor) # Data cleaning
# ==================== Statistical Modeling ====================
library(lme4) # Linear mixed-effects models
library(DHARMa) # Residual diagnostics
library(broom.mixed) # Tidy model outputs
# ==================== Visualization ====================
library(ggpubr) # Publication-ready plots
library(gt) # Publication-quality tables
cat("✓ Packages loaded successfully\n")## ✓ Packages loaded successfully
# Import buoyant weight data (skeletal mass measurements)
bw_raw <- read_csv(here("data", "fish_regen_buoyantweight (1).csv")) %>%
clean_names() %>%
# Remove outliers identified in QC
filter(coral_id != 45) %>%
filter(coral_id != 91)
# Import surface area data
sa <- read_csv(here("data", "fish_regen_mastersheet_wound closure_necrosis_sa - mastersheet.csv")) %>%
clean_names() %>%
mutate(coral_id = as.character(coral_id))
# Import tank metadata
tank_df <- read_csv(here("data", "Copy of fish_regen_mastersheet_wound closure_necrosis - mastersheet.csv")) %>%
clean_names() %>%
select(coral_id, tank) %>%
mutate(coral_id = as.character(coral_id))
cat("✓ Data imported successfully\n")## ✓ Data imported successfully
## - Buoyant weight data: 134 observations
## - Surface area data: 72 corals
## - Tank assignments: 72 corals
# Split data into initial and final measurements
bw_initial <- bw_raw %>%
filter(date == "initial") %>%
select(coral_id, wound, fish, initial_weight = bouyantweight_g)
bw_final <- bw_raw %>%
filter(date == "final") %>%
select(coral_id, final_weight = bouyantweight_g)
# Join into wide format (one row per coral)
bw_wide <- left_join(bw_initial, bw_final, by = "coral_id") %>%
mutate(
delta_mass = final_weight - initial_weight,
coral_id = as.character(coral_id)
)
cat("✓ Data reshaped to wide format\n")## ✓ Data reshaped to wide format
## - Initial weights: 67 corals
## - Final weights: 67 corals
# Adjust surface area to account for wound area removed
sa_adjusted <- sa %>%
select(-fish, -date_collected, -date_fragment) %>%
mutate(
# Subtract wound area: 3 cm² for large wounds, 1 cm² for small
sa_col = ifelse(wound == "Large", sa_cal - 3,
ifelse(wound == "Small", sa_cal - 1, sa_cal))
) %>%
select(-wound) # Remove to avoid duplication in join
# Merge buoyant weight with adjusted surface area
bw_sa_wide <- left_join(bw_wide, sa_adjusted, by = "coral_id") %>%
mutate(
wound = factor(wound, levels = c("No Wound", "Small", "Large")),
tank = as.factor(tank)
) %>%
mutate(tank = factor(tank, levels = sort(unique(as.numeric(as.character(tank))))))
cat("✓ Surface area data merged\n")## ✓ Surface area data merged
## - Mean initial SA: 69.38 cm²
Problem: Larger corals grow faster simply because they’re bigger, not necessarily because they’re healthier.
Solution: Allometric scaling adjusts growth for initial size using the relationship:
\[\log(\text{Final Mass}) = \alpha + b \times \log(\text{Initial Mass})\]
Where \(b\) is the allometric coefficient. We then calculate size-corrected growth as:
\[\text{Scaled Growth} = \log\left(\frac{\text{Final Mass}}{\text{Initial Mass}^b}\right)\]
# Create dataset with log-transformed weights
log_df <- bw_merged %>%
mutate(
log_i_weight = log(initial_weight),
log_f_weight = log(final_weight),
growth_ratio = final_weight / initial_weight
) %>%
drop_na(log_i_weight, log_f_weight, tank)
cat("✓ Allometric dataset prepared\n")## ✓ Allometric dataset prepared
## - Valid observations: 67 corals
# Linear mixed-effects model: final ~ initial + random(tank)
model_lmer <- lmer(log_f_weight ~ log_i_weight + (1 | tank), data = log_df)
cat("=== Allometric Model Summary ===\n")## === Allometric Model Summary ===
## Linear mixed model fit by REML ['lmerMod']
## Formula: log_f_weight ~ log_i_weight + (1 | tank)
## Data: log_df
##
## REML criterion at convergence: -515.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1853 -0.5898 0.0751 0.6736 2.0313
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 3.518e-07 0.0005931
## Residual 1.899e-05 0.0043574
## Number of obs: 67, groups: tank, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.044640 0.005354 8.338
## log_i_weight 0.990376 0.001961 504.931
##
## Correlation of Fixed Effects:
## (Intr)
## log_i_weght -0.994
# Extract allometric coefficient (b-value) with 95% CI
b_value <- tidy(model_lmer, effects = "fixed", conf.int = TRUE) %>%
filter(term == "log_i_weight")
cat("\n=== Allometric Coefficient (b-value) ===\n")##
## === Allometric Coefficient (b-value) ===
## # A tibble: 1 × 7
## effect term estimate std.error statistic conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed log_i_weight 0.990 0.00196 505. 0.987 0.994
# Generate predicted values
fit_df <- log_df %>%
select(log_i_weight) %>%
distinct() %>%
arrange(log_i_weight) %>%
mutate(predicted = predict(model_lmer, newdata = ., re.form = NA))
# Create allometric scaling plot
allometric_plot <- ggplot(log_df, aes(x = log_i_weight, y = log_f_weight)) +
geom_point(size = 2.2, shape = 21, fill = "steelblue", color = "black",
alpha = 0.8, stroke = 0.3) +
geom_line(data = fit_df, aes(x = log_i_weight, y = predicted),
color = "darkred", linewidth = 1.2) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed",
color = "gray60", linewidth = 0.8) +
labs(
title = "Allometric Scaling of Coral Buoyant Weight",
subtitle = "Mixed-effects model with tank as random intercept",
x = "Log Initial Weight (g)",
y = "Log Final Weight (g)"
) +
theme_pubr(base_size = 14)
print(allometric_plot)# Save figure
ggsave(
filename = here("figures", "growth", "allometric_scaling_plot.pdf"),
plot = allometric_plot,
width = 8, height = 6, dpi = 600, device = cairo_pdf
)
ggsave(
filename = here("figures", "growth", "allometric_scaling_plot.png"),
plot = allometric_plot,
width = 10, height = 7.5, dpi = 300, bg = "white"
)
cat("✓ Allometric plot saved\n")## ✓ Allometric plot saved
# Create formatted table of allometric coefficient
b_value_gt <- tidy(model_lmer, effects = "fixed", conf.int = TRUE) %>%
filter(term == "log_i_weight") %>%
transmute(
Term = "Allometric slope (b-value)",
Estimate = round(estimate, 3),
`95% CI` = paste0("[", round(conf.low, 3), ", ", round(conf.high, 3), "]")
)
b_value_gt %>%
gt() %>%
tab_header(
title = "Estimated Allometric Slope",
subtitle = "Mixed-effects model with tank as random intercept"
) %>%
cols_label(
Term = "",
Estimate = "Estimate",
`95% CI` = "95% Confidence Interval"
) %>%
tab_options(
table.font.size = 14,
heading.title.font.size = 16,
heading.subtitle.font.size = 13
)| Estimated Allometric Slope | ||
| Mixed-effects model with tank as random intercept | ||
| Estimate | 95% Confidence Interval | |
|---|---|---|
| Allometric slope (b-value) | 0.99 | [0.987, 0.994] |
# Create dataset with surface area and time adjustments
bw_sa_df <- bw_sa_wide %>%
drop_na(initial_weight, final_weight, wound) %>%
filter(coral_id != 103) %>% # Remove outlier
mutate(
delta_mass = final_weight - initial_weight,
wound = factor(wound, levels = c("No Wound", "Small", "Large"))
)
bw_sa_stats <- bw_sa_df %>%
drop_na(initial_weight, final_weight, fish, wound, tank) %>%
mutate(
wound = factor(wound, levels = c("No Wound", "Small", "Large")),
fish = factor(fish),
tank = factor(tank),
# Growth rate: mass change per surface area per day
bw_sa_time = (final_weight - initial_weight) / (sa_cal * 21)
)
cat("✓ SA-adjusted growth data prepared\n")## ✓ SA-adjusted growth data prepared
## - Valid observations: 66 corals
## - Mean growth rate: 0.00021 g/cm²/day
sa_growth_plot <- ggplot(bw_sa_stats, aes(x = sa_cal, y = delta_mass, color = wound)) +
geom_point(size = 2.2, alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE, linewidth = 1.1) +
scale_color_brewer(palette = "Dark2") +
labs(
title = "Coral Growth as a Function of Surface Area",
subtitle = "Colored by Wound Treatment",
x = "Surface Area (cm²)",
y = "Mass Change (g)",
color = "Wound Treatment"
) +
theme_pubr(base_size = 14)
print(sa_growth_plot)# Save figure
ggsave(
filename = here("figures", "growth", "sa_growth_relationship.pdf"),
plot = sa_growth_plot,
width = 8, height = 6, dpi = 600, device = cairo_pdf
)
ggsave(
filename = here("figures", "growth", "sa_growth_relationship.png"),
plot = sa_growth_plot,
width = 10, height = 7.5, dpi = 300, bg = "white"
)
cat("✓ SA-growth plot saved\n")## ✓ SA-growth plot saved
# Extract b-value from allometric model
b_est <- b_value$estimate[1]
# Calculate allometrically-scaled growth
growth_stats <- bw_merged %>%
drop_na(initial_weight, final_weight, fish, wound, tank) %>%
mutate(
wound = factor(wound, levels = c("No Wound", "Small", "Large")),
fish = factor(fish),
tank = factor(tank),
# Size-corrected growth using estimated b-value
growth_scaled = log(final_weight / initial_weight^b_est)
)
cat("✓ Scaled growth variable calculated\n")## ✓ Scaled growth variable calculated
## - Using b-value: 0.99
## - Mean scaled growth: 0.0446
# Full model with interaction
mod_int <- lmer(growth_scaled ~ fish * wound + (1 | tank), data = growth_stats)
# Additive model (no interaction)
mod_add <- lmer(growth_scaled ~ fish + wound + (1 | tank), data = growth_stats)
# Fish-only model
mod_fish <- lmer(growth_scaled ~ fish + (1 | tank), data = growth_stats)
# Wound-only model
mod_wound <- lmer(growth_scaled ~ wound + (1 | tank), data = growth_stats)
# Null model
mod_null <- lmer(growth_scaled ~ 1 + (1 | tank), data = growth_stats)
cat("✓ All models fitted successfully\n")## ✓ All models fitted successfully
# Test for interaction
lrt_interaction <- anova(mod_add, mod_int, test = "Chisq")
# Test for fish effect
lrt_fish <- anova(mod_wound, mod_add, test = "Chisq")
# Test for wound effect
lrt_wound <- anova(mod_fish, mod_add, test = "Chisq")
# Compile results
lrt_table <- tibble(
Test = c("Interaction (Fish × Wound)", "Fish Effect", "Wound Effect"),
ChiSq = c(lrt_interaction$Chisq[2], lrt_fish$Chisq[2], lrt_wound$Chisq[2]),
Df = c(lrt_interaction$Df[2], lrt_fish$Df[2], lrt_wound$Df[2]),
p_value = c(lrt_interaction$`Pr(>Chisq)`[2], lrt_fish$`Pr(>Chisq)`[2], lrt_wound$`Pr(>Chisq)`[2])
) %>%
mutate(across(where(is.numeric), round, 3))
# Format table
lrt_table %>%
gt() %>%
tab_header(
title = "Likelihood Ratio Tests: Allometrically-Scaled Growth"
) %>%
cols_label(
Test = "Model Comparison",
ChiSq = "χ²",
Df = "df",
p_value = "P-value"
) %>%
fmt_number(columns = c(ChiSq, p_value), decimals = 3) %>%
tab_options(
table.font.size = 14,
heading.title.font.size = 16
)| Likelihood Ratio Tests: Allometrically-Scaled Growth | |||
| Model Comparison | χ² | df | P-value |
|---|---|---|---|
| Interaction (Fish × Wound) | 4.225 | 2 | 0.121 |
| Fish Effect | 0.011 | 1 | 0.918 |
| Wound Effect | 1.856 | 2 | 0.395 |
# Full model with interaction
mod_int_sa <- lmer(bw_sa_time ~ fish * wound + (1 | tank), data = bw_sa_stats)
# Additive model
mod_add_sa <- lmer(bw_sa_time ~ fish + wound + (1 | tank), data = bw_sa_stats)
# Fish-only
mod_fish_sa <- lmer(bw_sa_time ~ fish + (1 | tank), data = bw_sa_stats)
# Wound-only
mod_wound_sa <- lmer(bw_sa_time ~ wound + (1 | tank), data = bw_sa_stats)
# Null
mod_null_sa <- lmer(bw_sa_time ~ 1 + (1 | tank), data = bw_sa_stats)
cat("✓ Growth rate models fitted\n")## ✓ Growth rate models fitted
# Likelihood ratio tests
lrt_interaction_sa <- anova(mod_add_sa, mod_int_sa, test = "Chisq")
lrt_fish_sa <- anova(mod_wound_sa, mod_add_sa, test = "Chisq")
lrt_wound_sa <- anova(mod_fish_sa, mod_add_sa, test = "Chisq")
# Compile results
lrt_table_sa <- tibble(
Test = c("Interaction (Fish × Wound)", "Fish Effect", "Wound Effect"),
ChiSq = c(lrt_interaction_sa$Chisq[2], lrt_fish_sa$Chisq[2], lrt_wound_sa$Chisq[2]),
Df = c(lrt_interaction_sa$Df[2], lrt_fish_sa$Df[2], lrt_wound_sa$Df[2]),
p_value = c(lrt_interaction_sa$`Pr(>Chisq)`[2], lrt_fish_sa$`Pr(>Chisq)`[2], lrt_wound_sa$`Pr(>Chisq)`[2])
) %>%
mutate(across(where(is.numeric), round, 3))
# Format table
lrt_table_sa %>%
gt() %>%
tab_header(
title = "Likelihood Ratio Tests: Growth Rate (g/cm²/day)"
) %>%
cols_label(
Test = "Model Comparison",
ChiSq = "χ²",
Df = "df",
p_value = "P-value"
) %>%
fmt_number(columns = c(ChiSq, p_value), decimals = 3) %>%
tab_options(
table.font.size = 14,
heading.title.font.size = 16
)| Likelihood Ratio Tests: Growth Rate (g/cm²/day) | |||
| Model Comparison | χ² | df | P-value |
|---|---|---|---|
| Interaction (Fish × Wound) | 3.051 | 2 | 0.218 |
| Fish Effect | 0.058 | 1 | 0.809 |
| Wound Effect | 0.228 | 2 | 0.892 |
# Prepare data with correct factor ordering
bw_sa_stats$wound <- factor(bw_sa_stats$wound, levels = c("No Wound", "Small", "Large"))
bw_sa_stats$fish <- factor(bw_sa_stats$fish, levels = c("No Fish", "Fish"))
# Create publication-quality figure
growth_rate_plot <- bw_sa_stats %>%
ggplot(aes(x = wound, y = bw_sa_time, color = fish, shape = fish)) +
# Individual data points
geom_jitter(
aes(color = fish),
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.6),
size = 2.5,
alpha = 0.25,
shape = 16
) +
# Mean points
stat_summary(
fun = mean,
geom = "point",
size = 5,
position = position_dodge(width = 0.6)
) +
# Error bars (mean ± SD)
stat_summary(
fun.data = function(x) {
data.frame(y = mean(x), ymin = mean(x) - sd(x), ymax = mean(x) + sd(x))
},
geom = "errorbar",
width = 0,
size = 1,
position = position_dodge(width = 0.6)
) +
# Scales
scale_color_manual(
values = fish_palette,
name = "Fish Presence",
labels = c("No Fish", "Fish")
) +
scale_shape_manual(
values = c(16, 18),
name = "Fish Presence",
labels = c("No Fish", "Fish")
) +
# Theme
theme_minimal(base_size = 16) +
theme(
panel.grid = element_blank(),
axis.title = element_text(size = 20, face = "bold"),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 16),
legend.position = "top",
legend.title = element_text(size = 18, face = "bold"),
legend.text = element_text(size = 16),
legend.key.size = unit(1.5, "lines"),
legend.margin = margin(0, 0, 10, 0),
legend.box.spacing = unit(0.5, "lines"),
axis.line.x = element_line(color = "black", size = 0.5),
axis.line.y = element_line(color = "black", size = 0.5)
) +
labs(
y = "Growth Rate (mg/cm²/day)",
x = "Wound Size"
)
print(growth_rate_plot)# Save figure
ggsave(
filename = here("figures", "growth", "summary_line_growth_rate_by_wound_fish.pdf"),
plot = growth_rate_plot,
width = 7, height = 5, dpi = 600, device = cairo_pdf
)
ggsave(
filename = here("figures", "growth", "summary_line_growth_rate_by_wound_fish.png"),
plot = growth_rate_plot,
width = 8, height = 6, dpi = 300, bg = "white"
)
cat("✓ Growth rate figure saved\n")## ✓ Growth rate figure saved
## === Diagnostics for Growth Rate Model (bw_sa_time) ===
sim_growth_full <- simulateResiduals(fittedModel = mod_int_sa, n = 1000)
# Save diagnostic plot
png(
here("figures", "diagnostics", "growth_dharma_full.png"),
width = 7, height = 6, units = "in", res = 600, bg = "white"
)
plot(sim_growth_full, main = "DHARMa Diagnostics: Growth Rate Full Model")
dev.off()## quartz_off_screen
## 2
##
## Dispersion test:
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91737, p-value = 0.67
## alternative hypothesis: two.sided
##
## Outlier test:
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_growth_full
## outliers at both margin(s) = 0, observations = 66, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.00000000 0.05435885
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0
bw_sa_stats$pred_growth <- predict(mod_int_sa)
p_growth_fit <- ggplot(bw_sa_stats, aes(x = pred_growth, y = bw_sa_time)) +
geom_point(aes(color = fish), size = 2.5, alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linewidth = 1) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
scale_color_manual(values = fish_palette) +
labs(
x = "Predicted Growth Rate (mg/cm²/day)",
y = "Observed Growth Rate (mg/cm²/day)",
title = "Growth Model Fit: Predicted vs Observed",
color = "Fish Presence"
) +
theme_bw(base_size = 14) +
theme(legend.position = "top")
print(p_growth_fit)##
##
## === Diagnostics for Allometric Scaling Model ===
sim_allometric <- simulateResiduals(fittedModel = model_lmer, n = 1000)
png(
here("figures", "diagnostics", "allometric_dharma.png"),
width = 7, height = 6, units = "in", res = 600, bg = "white"
)
plot(sim_allometric, main = "DHARMa Diagnostics: Allometric Scaling Model")
dev.off()## quartz_off_screen
## 2
##
## Dispersion test:
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98124, p-value = 0.972
## alternative hypothesis: two.sided
##
## Outlier test:
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_allometric
## outliers at both margin(s) = 1, observations = 67, p-value = 0.1254
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0003778063 0.0803764506
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01492537
##
## === Growth Model Diagnostics Summary ===
## Residual standard deviation: 6.533866e-05
## Number of observations: 66
##
## Residual summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.674e-04 -3.278e-05 -1.770e-06 0.000e+00 3.712e-05 1.963e-04
## R version 4.5.1 (2025-06-13)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.0.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
##
## time zone: America/Los_Angeles
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gt_1.0.0 ggpubr_0.6.1 broom.mixed_0.2.9.6
## [4] DHARMa_0.4.7 lme4_1.1-37 Matrix_1.7-3
## [7] janitor_2.2.1 here_1.0.1 lubridate_1.9.4
## [10] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
## [13] purrr_1.1.0 readr_2.1.5 tidyr_1.3.1
## [16] tibble_3.3.0 ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.4 rlang_1.1.6 magrittr_2.0.3 snakecase_0.11.1
## [5] furrr_0.3.1 compiler_4.5.1 mgcv_1.9-3 systemfonts_1.2.3
## [9] vctrs_0.6.5 pkgconfig_2.0.3 crayon_1.5.3 fastmap_1.2.0
## [13] backports_1.5.0 labeling_0.4.3 utf8_1.2.6 promises_1.3.3
## [17] rmarkdown_2.29 tzdb_0.5.0 nloptr_2.2.1 ragg_1.4.0
## [21] bit_4.6.0 xfun_0.53 cachem_1.1.0 jsonlite_2.0.0
## [25] later_1.4.3 broom_1.0.9 parallel_4.5.1 R6_2.6.1
## [29] gap.datasets_0.0.6 bslib_0.9.0 stringi_1.8.7 qgam_2.0.0
## [33] RColorBrewer_1.1-3 parallelly_1.45.1 car_3.1-3 boot_1.3-31
## [37] jquerylib_0.1.4 Rcpp_1.1.0 iterators_1.0.14 knitr_1.50
## [41] httpuv_1.6.16 splines_4.5.1 timechange_0.3.0 tidyselect_1.2.1
## [45] abind_1.4-8 yaml_2.3.10 doParallel_1.0.17 codetools_0.2-20
## [49] listenv_0.9.1 lattice_0.22-7 plyr_1.8.9 shiny_1.11.1
## [53] withr_3.0.2 evaluate_1.0.4 future_1.67.0 xml2_1.4.0
## [57] pillar_1.11.0 gap_1.6 carData_3.0-5 foreach_1.5.2
## [61] reformulas_0.4.1 generics_0.1.4 vroom_1.6.5 rprojroot_2.1.0
## [65] hms_1.1.3 scales_1.4.0 minqa_1.2.8 globals_0.18.0
## [69] xtable_1.8-4 glue_1.8.0 tools_4.5.1 ggsignif_0.6.4
## [73] grid_4.5.1 rbibutils_2.3 nlme_3.1-168 Formula_1.2-5
## [77] cli_3.6.5 textshaping_1.0.1 gtable_0.3.6 rstatix_0.7.2
## [81] sass_0.4.10 digest_0.6.37 farver_2.1.2 htmltools_0.5.8.1
## [85] lifecycle_1.0.4 mime_0.13 bit64_4.6.0-1 MASS_7.3-65
Analysis Complete ✓
All figures and diagnostics saved to figures/
directory.